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We calculate the single-particle spectral function for doped bilayer graphene in the low energy 
limit, described by two parabolic bands with zero band gap and long range Coulomb interaction. 
Calculations are done using thermal Green's functions in both the random phase approximation 
(RPA) and the fully self-consistent GW approximation. RPA (in line with previous studies) yields 
a spectral function which apart from the Landau quasiparticle peaks shows additional coherent fea- 
tures interpreted as plasmarons, i.e. composite electron-plasmon excitations. In GW the plasmaron 
becomes incoherent and peaks are replaced by much broader features. The deviation of the quasipar- 
ticle weight and mass renormalization from their non-interacting values is small which indicates that 
bilayer graphene is a weakly interacting system. The electron energy loss function, Im[— e~ (u)] 
shows a sharp plasmon mode in RPA which in GW approximation becomes less coherent and thus 
consistent with the weaker plasmaron features in the corresponding single-particle spectral function. 

PACS numbers: 



I. INTRODUCTION 

Since its fabrication, graphene [IH5] (a single layer of 
graphite) has been of interest for both theoreticians and 
experimentalists. It is a two-dimensional (2D) crystal 
with carbon atoms arranged on a honeycomb lattice with 
two sublattices. Due to its unique properties(e.g. high 
mobility even in highly doped cases) it opens new per- 
spectives for engineering and is a candidate material for 
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FIG. 1: (Color online) Single particle spectral function for 
bilayer graphene in the low energy limit at r s = 3. (a) - RPA, 
(b) - GW. The bare bands eg = ±fc 2 (in units k F — 1, e F — 1) 
are rotationally symmetric (the patchy appearance is due to 
the finite k-space resolution), (c) and (d) are the cuts (dash- 
dotted lines) in (a) and (b), respectively. Dashed lines are 
guides to the eye for plasmaron dispersions. 



future nanoelectronic and spintronic devices [5]. 

The subject of this article is the closely related bi- 
layer graphene, formed by stacking two graphene layers 
in Bernal " AB" stacking sequence in which the two layers 
are rotated by 60 degrees. These are coupled by inter- 
layer tunneling between A and B sublattice sites with 
the hopping parameter t± m 0.39 eV [7j. 

Bilayer graphene shares some features with both 
graphene and the ordinary two-dimensional electronic 
gas (2DEG). Its dispersion is quadratic, similar to a 
2DEG but the effective Hamiltonian is chiral with zero 
band gap as in the case of graphene [7J [5] • In both sin- 
gle layer and bilayer graphene the charge carrier density 
can be controlled by application of a gate voltage, a fun- 
damental effect for potential technological applications 
[Hinj-In addition, for bilayer graphene even the band gap 
is tunable with great potential for device applications 

pun]. 

Apart from the dispersion relation, the property which 
makes bilayer graphene different from that of a single 
layer is its coupling parameter being a function of the 
carrier density r s ~ n -1 / 2 [21 QT]. In other words the 
strength of Coulomb interaction is tunable, while the cou- 
pling parameter for the single layer graphene is constant 
r s ~ n° and lies in the interval < r s < 2.2. By com- 
paring the values of r s for single- and bilayer graphene 
(r s ~ 68.5 x 10 5 / \/ri, where n is the number of carriers 
per cm with n m 10 9 — 5 x 10 12 cm~ 2 ) in vacuum it 
is clear that the strength of the Coulomb interaction can 
be much larger in bilayer graphene [2]. 

The electronic structure of bilayer graphene [Jj [12] 
is characterised by the single particle spectral function 
AAui), which can be measured experimentally by angle 
resolved photo-emission spectroscopy (ARPES)[T31 [T3] . 
It obeys the sum rule J ^AAu)) = 1 and can be in- 
terpreted as the probability distribution of an electron 
having momentum k and energy uj. Sensarma et al. |llj 
studied how Coulomb interaction affects the single parti- 
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cle spectral function of bilayer graphene away from half- 
filling. The authors used RPA to calculate that doped 
bilayer graphene is a Fermi liquid in the low energy limit, 
with a sharp quasiparticle peak. They also found addi- 
tional weaker peak structures that they interpreted as 
plasmarons; a quasiparticle formed by the coupling be- 
tween electron and plasmon, as originally predicted by B. 
Lundqvist [15] , Studying the physics of interaction be- 
tween electrons and plasmons in graphene is particularly 
interesting because of recently proposed " plasmonic" de- 
vices that could merge photonics and electronics [14] . 

Experimentally plasmarons in the single layer 
graphene were observed by A. Bostwick et al [14] us- 
ing angle-resolved photoemission spectroscopy. Apart 
from the two single particle crossing bands, two addi- 
tional bands were observed and interpreted as a spec- 
trum of plasmarons. The experimentally measured spec- 
tral function compares qualitatively with that obtained 
within RPA. 

In this paper we compute numerically the single- 
particle spectral function A?(uS) for doped bilayer 
graphene in the low energy two-band approximation in 
both RPA and the fully self-consistent GW approxima- 
tion [TBHH]. We use a thermal Green's function formal- 
ism, based on a finite set of imaginary frequencies and 
analytic continuation to real frequencies for the single- 
particle Green's functions in a controlled manner [19 . 

The results, presented in Fig. [lj show the spectral 
function with long lived Landau quasiparticles and satel- 
lite plasmaron peaks in RPA (Fig. fl| (a)) and confirm 
the results of analytic calculations [TUTTl |2"U] , whereas in 
the GW approximation the plasmaron peaks are replaced 
by broad shoulders (Fig. [I] (b)). It has been emphati- 
cally argued that self-consistent GW underestimates the 
coherence of collective excitations [HI l2TT - f23] and our 
results showing a marked difference between the satel- 
lite peaks in RPA and GW most likely agree with this. 
Nevertheless we argue that the GW results are valuable 
as a benchmark for more sophisticated self-consistent ap- 
proaches including vertex corrections to the polarization. 

Below we describe our calculations in more detail. 

II. GW APPROXIMATION 

The GW approximation is derived perturbatively from 
the Hedin's equations jTHJ [IT] , giving the self-energy 
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where W^{iuj n ) and G^{iuj n ) are the dressed interaction 
and Green's function, respectively (all quantum numbers, 
such as momentum, spin, etc., are incorporated in k and 
q). The argument of the Green's function is the fermionic 
Matsubara frequency uj n — (2n + l)ir/f3, whereas the 



dressed interaction is a function of the bosonic Matsub- 
ara frequency ui n = 2mr/f3 with n integer. After com- 
puting T,'3 w (iuj n ) (first diagram in Fig. [5| one should, 
in general, add the Hartree diagram (second diagram in 
Fig. [2| to it which in case of long-range Coulomb inter- 
action gives zero contribution because it is cancelled by 
the positive background charge |24] , The approximation 




FIG. 2: Self-energy in the GW approximation. Double wiggly 
line, single wiggly line and double line correspond to dressed 
interaction, bare interaction and dressed Green's function, re- 
spectively 

has the same form as the standard Hartree-Fock (HF) 
approximation. The difference is that the latter uses the 
bare Green's function and interaction, while the former 
is based on the dressed Green's function Gj:(uj) and dy- 
namically screened (dressed) interaction. 

The screened interaction W^{iu n ) is an infinite geo- 
metric series of diagrams (Fig. [3| consisting of the bare 
interaction V q and the irreducible polarization diagram 
H$(iuj n ) which in GW is given by 

\ ' m— — oo 

(2) 

where g is the degeneracy factor. In RPA it is com- 
puted using the same relation but with the full Green's 
functions replaced by the bare ones. So, the RPA polar- 
ization is just the zeroth order term in the expansion of 
ilg-(w„) in the bare interaction. After summing up the 
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FIG. 3: Screened interaction in the GW approximation is 
given by geometric series. Bubble diagram represents polar- 
ization Yl^(iu)n). 

geometric series one obtains the following expression for 
the screened interaction 
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The effective bare Coulomb interaction for bilayer 

9 2 

graphene is given by V q = where k represents the 
background dielectric constant. [12] Using Ep and kp as 
units of energy and momentum, respectively enables us 
to write V q in terms of the dimensionless coupling pa- 
rameter r s = e 2 gm/{kpK): 

nr. 



V„ 



(4) 



In the non- interacting limit T\^{iw n ) can be written in 
the following simple form: 



U°Mu n ) = -g 



E 

s.s' 



d 2 k (f s k -fl' +q )F s ,Ak,k + g) 
(2tt) 2 iuj 



i 1 



(10) 

where ft = 1/(1 + e^ ek ^) is the Fermi distribution 
function and 



III. BILAYER GRAPHENE 



F,,./(£, k + q) = -Jr(l + sa % )(l + s'a %+q ) 



A. Effective model 



-(l + SS 'cos(2^ + ,)) 



The low energy limit of bilayer graphene is valid if 
the scale of all relevant energies are smaller than the in- 
terlayer hopping parameter t± such that the two outer 
bands can be ignored. Incorporating the two layers as 
an additional index, the Hamiltonian can be represented 
as a 2x2 matrix and thus four-band model is reduced to 
the effective two-band model with the total degeneracy 
of g = 4 (due to the spin and valley index) , 

1_ f (k x 

2m \ (kx - iky) 2 

It is clear that the corresponding energy spectrum is 
parabolic, 

k 2 

(6) 
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£fc = ± 



2m 



m = t±/(2v F ) « 0.054m e is the effective mass of the elec- 
tron in the low energy limit, m e being the free electron 
mass. Corresponding free Fermionic Matsubara Green's 
function is 



G°JiLLi n ) = (iw r , 



2 ' iu). 

s=± 



s|e fc | + /i' 
(7) 

s indexes the conductance and valence bands, fi is the 
chemical potential and <7g is given by 



k 2 
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IL q -(iw n ) = -- / -—- J2 Tr(Gj:(iu m )x 
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xG k+q( iuJ '"< + iuJ rn))- (9) 



with Oj: g , - being the angle between the vectors k and k+ 
q. Tlq-(iuj n ) is an angle-independent function which can be 
seen by extracting the angle 9q- using the rotation of the 
integration variable in Eq. 9] with 9^. Consequently, the 
screened interaction is angle-independent as well. Note 
that since the polarization is a scalar due to the trace 
in Eq. [9] the screened interaction remains to be a scalar 
quantity as well. 

Clearly, the GW self-energy is generalised to 



S ? (iw n ) 



d 2 q 
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„ (11) 

In order to the see the matrix structure of the T,^ w (iu> n ) 
we perform t he following integration variable transforma- 
tions in Eq. jnjwith replaced by GS _, 



and 



q x =k-q 



<f 2 = TZ(ir + %)<fi, 



(12) 



(13) 



where 1Z(tt + 6%) denotes the rotation matrix with the 
angle tt + dz- Therefore the self-energy can be rewritten 

as 



where k± = k x ± iky, a± — (a± ± i&2)/2, o\ and a-x with 
are Pauli matrices, 9^ is the angle of the vector k with 
respect to the x-axis. 

Let us rewrite Eq. [2] for the two-band model. After 
summing over all internal indicies (which in this case is anc j 
just the matrix index of G^(iw n ) ) when computing the 
polarization diagram we obtain the following expression 
for II ? (iw„), 



E = 



1 



S S+e i29 £ 



d 2 q 
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(14) 



(15) 



(16) 



Now, if one makes the variable transformation 9^ 2 = — % 2 
in S + (£_) it becomes obvious that E + = £_ which 
means that fl9 (iui n ) and the fully interacting Green's 

function GAiu) n ) have and retain the same structure as 
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the one of the free Green's function GS(iu;„) (Eq. |l7j|l8| 
throughout the whole self-consistent calculation. 



periodic set of the Green's function values in the Mat- 
subara frequency space 



G%[ioj n ) = 



(17) 



G° k (iuj n ) = r]Cothri(iu! n - e k ). 



(20) 



af = (iu n - \e k \ ± (w n + \e k \ + fj,)' 



So, it is sufficient to set up calculations for £^ (iw„) 

and G^(iu! n ) only at 9g = 0. 

We start the self-consistent calculation by discretising 
momenta and angles. Since our interest is focused on 
the low energy properties the absolute value of k ranging 
from to 4 is discretised into 40 points logarithmically 
giving denser number of points around kp. The rest of 



Here r\ = and e k is a single-particle excitation spec- 

(18) 

trum for a given model. It is obvious that Eq. [20] is 
periodic under iu n — > iuj n + ii^Ni = —< The peri- 



odized full Green's function is given by 

G k (iuj n ) = ?7C0thr?(ia;„ - e k - T, k (ioj n )), 



(21) 



which has the correct non-interacting limit and together 
with G k (iui n ) yields standard continuum expression for 
the Greens function as N tends to oo ( rj — > ). Due to 



9j£ and 



the integration variables ( 
linearly. \q\ is discretized into 80 points and lies in the 
interval [1/80, 4] while the number of the discretization 
points for and 0$ is 10. First, the free Greens's function 
is evaluated at = and then it is rotated by an angle 
9j: in order to obtain the polarization (Eq. Wm. Then the 
screened interaction is computed using Eq.TS which en- 
ables us to evaluate the GW self-energy (Eq. II]). After 
calculating Y,9 W (iuj n ) at 9^ = we update the Green's 
function through Eq. [22j This is done repeatedly: if the 
procedure converges to a fixed point, a solution has been 
found. The calculations are done at T/ep = 1/10 with 
N = 121 number of Matsubara frequencies. 



B. Periodized Green's functions 

The GW approximation implies a self-consistent nu- 
merical calculation, which may be solved iteratively |27j . 
Obviously these calculations include very demanding op- 
erations including infinite sums over Matsubara frequen- 
cies. In order to cope in numerical calculations with these 
kinds of problems we use a new formalism for finite tem- 
perature fermionic thermal Green's functions in the sin- 
gle band case described in [THj and summarised below. 

Performing numerical calculations using thermal 
Green's functions [25J [2H| may be done by the discretiza- 
tion of imaginary time. Since the fermionic thermal 
Green's function is anti-periodic over r G [— f3, f3] domain 
with the period /? we discretize the interval t £ [0, f3] into 
N evenly spaced points, r = ((3/N)j 7 j — 1, ...,N — 1. 
Due to the discontinuity of the fermionic Green's func- 
tion at t = (limits r — > 0~ and r — > + differ 
from each other) some specific value must be assigned 
to G k (rj = 0) when doing numerical computations. We 
define G k {r.j = 0) by the average of G k {r — _ ) and 
Gfe(r = + ). After applying discrete Fourier transforma- 
tion to the non-interacting thermal Green's function 

G° k (r) = e-^[{n k - l)9(r) + n k 9(-r)} (19) 

where n k = (ct c k ) is the occupation number we obtain 



-) are discretized ^ e non t r i v i a l hyperbolic function in Eq. 21 one can not 



define the self-energy using simply Gq 1 and G _1 as it is 
done in the standard theory. In this case the self-energy 
is defined by the amputated skeleton diagrams ([5SJ, see 



Sec. 5.1) and through Eq. 21 



In the case of two band model G k (iuj n ) is generalised 



to 



G s (iuj n ) =r i cothr,{{Gi(ioj n ))- 1 - ±^ n )), (22) 



where 



= (iu> n + y)¥ - \tk\o k - 



The periodized Green's function for both single and 
two-band cases is consistent with the corresponding 
Luttinger-Ward T-functional [531 [3D] (the former is con- 
sistent with the T-functional as presented in Eq. 4 in [HI] 
while the former - with the same equation where G k (iu> n ) 
is replaced by G^{iiu n )). 

To perform analytic continuation for G k (ioj n ) we first 
rewrite it by means a conformal transformation in a new 
basis where it can be represented as a sum of simple 
poles. Then the Pade method [3T] of fitting to a rational 
function is used which enables us to evaluate the Green's 
function on the real frequency axis. In the case of a two- 
band model the trace of G^(iu: n ) is used as an input to 
the same procedure of analytic continuation as the one 
carried out for G k (iu> n ). 



C. Spectral function 



The spectral function is given by 
1 



Im[TrGAu + i0 + )}. 



(23) 



where we perform analytic continuation after applying 
trace to Gg(iw„). In order to study the low energy prop- 
erties we also compute the spectral function projected on 
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the conductance band 



A % {s = +,u) = —Im[G £ (s = +,uj + i0+)], (24) 



where Gj:(s = represents the eigenvalue of Gg(iu) n ) 

corresponding to the upper band after analytic continu- 
ation to the real axis. In Fig. [I] and [5] we present the 
spectral functions (left column) for different values of kp 
together with the corresponding self-energies (right col- 
umn) in the GW approximation and RPA at r s — 3 and 
r s = 7, respectively. 

As the plots show the spectral weight in the RPA away 
from k F has two peaks: the main Landau quasiparticle 
peak and plasmaron peaks. The presence of the plas- 
maron excitation also give jumps in the real and imagi- 
nary parts of the corresponding self-energies. The RPA 
plasmaron excitation has lower weight at r s = 7 than the 
one at r s = 3, although the spectral functions have qual- 
itatively same behaviour which is also noticeable in the 
case of the GW approximation. Most of the structure 
obtained in the RPA is not presented in the GW ap- 
proximation. We interpret this as being due to stronger 
screening in GW. 

In Fig. [6]Ja) the electron energy loss spectrum 
Im[— e^ 1 (w)] (e q (ui) = 1 + V q H^(uj) - dielectric function) 
in RPA is plotted showing the plasmaron dispersion re- 
lation (black color) which is in a quite good agreement 
for small g-values with its analytic version (solid line) 
expanded up to second order in q [21 [20] , 



gE F q 



1 - 



r s q 

8k f 



Im[— e~ l (u))] was also calculated in the GW approxima- 
tion (Fig. [6jb) ) where the plasmon mode is less coherent 
than that in RPA which is in agreement with the fact 
that the plasmaron features in the GW spectral function 
are weaker than in RPA. 



D. Quasiparticle weight and effective mass 

The quasiparticle weight Z and renormalized mas m* , 
given in Table I, are computed for both the GW and RPA 
approximations using the formulas: 



Z = 

# 

m 
m 



1 



z- 1 



m dReY, k {u)=e F ) i 
k F dk 



1 



(25) 



(26) 



As expected, the quasiparticle weight decreases with 
increasing interaction strength because the interaction 
shifts the weight from the coherent quasiparticle peak 
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FIG. 4: (Color online) r s — 3. Left column: spectral weight 
in RPA (dashed line) and GW approximation (solid line) at 
k w 0.76fc F (a), k = k F (b), k « 1.20fc F (c). Right column: 
the real (blue solid line) and imaginary (red dashed line) part 
of the self-energy at k ~ 0.76fcF (d), k — kp (e), k w 1.20/cf 
(f) in RPA (thin line) and GW approximation (thick line). 
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FIG. 5: (Color online) r 3 — 7. Left column: spectral weight 
in RPA (dashed line) and GW approximation (solid line) at 
k ~ 0.76kp (a), k = kp (b), k ~ 1.20/cf (c). Right column: 
the real (blue solid line) and imaginary (red dashed line) part 
of the self-energy at k ~ 0.76fcF (d), k — kp (e), k w 1.20/cf 
(f) in RPA (thin line) and GW approximation (thick line). 



through incoherent scattering. Since the GW approxi- 
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FIG. 6: (Color online) /mf-e^ 1 (u))) in RPA (a) and GW ap- 
proximation (b) at r s — 7 (same color intensity scale on both 
plots). Green solid line in (a) represents the plasmon disper- 
sion expanded up to the second order in q. The unexpected 
discontinuities are artificial and due to difficulties with the 
analytic continuation of a two-particle function. 



mation does not yield the plasmaron peaks and the in- 
teraction gets more screened, most of the weight is con- 
centrated in the Landau quasiparticle which results in a 
bigger quasiparticle weight than that in the case of RPA. 
The mass renormalization is less than 7% in both ap- 
proximations meaning that we are dealing with a weakly 
interacting system. By comparing our results with the 



RPA 0.798 0.978 
GW 0.851 0.946 



Z m* /m 
RPA 0.685 0.986 
GW 0.806 0.929 



TABLE I: Quasiparticle weight Z and effective mass relative 
to the one of the free electron m*/m at r s = 3 (left) and 
r s = 7 (right). 

ones presented in [TT] one can see that the agreement is 
quite good. 



IV. CONCLUSION 

We present the single particle spectral function and 
self-energy for bilayer graphene in the low energy limit as 
described with a two band model. Calculations are done 
in both RPA and self-consistent GW using a discretized 
thermal Green's function formalism. In RPA, the spec- 
tral function and energy loss spectrum show prominent 
plasmaron peaks and sharp plasmon mode, respectively 
whereas in GW the plasmaron peaks are replaced by 
broad shoulders which is consistent with a less coherent 
plasmon mode. The RPA spectral function, quasiparticle 
weight and effective mass are in a good agreement with 
those in [11] computed using a conventional Matsubara 
Green's function method. 
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